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ABSTRACT 

One of the most puzzling facts about cuprate high-temperature superconductors in the lightly doped regime is the coexistence 
of uniform superconductivity and/or antiferromagnetism with many low-energy charge-ordered states in a unidirectional charge 
density wave or a bidirectional checkerboard structure. Recent experiments have discovered that these charge density waves 
exhibit different symmetries in their intra-unit-cell form factors for different cuprate families. Using a renormalized mean-field 
theory for a well-known, strongly correlated model of cuprates, we obtain a number of charge-ordered states with nearly 
degenerate energies without invoking special features of the Fermi surface. All of these self-consistent solutions have a pair 
density wave intertwined with a charge density wave and sometimes a spin density wave. Most of these states vanish in the 
underdoped regime, except for one with a large d-form factor that vanishes at approximately 19% doping of the holes, as 
reported by experiments. Furthermore, these states ceuld be modified te have a global superconducting order, with a nodal-like 
density of states at low energy. 


Introduction 

Ever since the discovery of the high-Tc superconductivity, many low-energy charge-ordered states in the cuprate have been 
discovered. Neutron scattering experimentsfl] first emphasised the doping dependence of incommensurate magnetic peaks 
associated with unidirectional magnetic patterns or stripes. Later, soft X-ray scattering[2] also confirmed the presence of charge 
orders with these stripes. However, these experiments were performed on the 2 \ 4 (La 2 -xSrxCu 04 ) cuprate family. For other 
cuprate families, the evidence for bond-centred unidirectional domains was found via scanning tunneling spectroscopy[3,4]. 
The charge density wave(CDW) order was also found to be induced by the external magnetic field[5]. Recently, more results 
regrading charge-ordered states [6-10], and electron-doped cuprates[l 1] have been reported. The periods of these CDW and 
their doping dependence are quite different for different cuprate families [10]. In addition to the unidirectional stripe pattern, 
some experiments have also reported the possible existence of a bidirectional charge-ordered checkerboard pattern[12,13]. The 
unidirectional charge-ordered states or stripes were found to have a dominant d-like symmetry for the intra-unit-cell form factor, 
measured on the two oxygen sites by using the resonant elastic x-ray scattering method[14,15] and via scanning tunneling 
spectroscopy (STS)[16]. However, different families seem to prefer different symmetries[14,15]. In the STS experiments[17], 
the density waves disappeared above 19% hole doping. Furthermore, the observation of these CDW states having nodal-like 
local density of states(LDOS) at low energy but strong spatial variation at high energy in STS[3] strongly implies a new 
unconventional superconducting state. 

The existence of these great varieties of charge-ordered states has created a great debate regarding whether the strong 
coupling Hubbard model or the t — J model[ 18] is the proper basic Hamiltonian to describe the cuprates. Many believe that 
these states “compete” with the superconductivity[19] and that their origin may reveal the fundamental understanding of the 
mechanism of high superconducting temperatures in cuprates. The recent detection of the d-form factor at an oxygen site 
instead of at a Cu site[14-16] also raises the question about the suitability of the effective one-band Hubbard or t — J model and 
the validity of replacing the oxygen hole with a Zhang-Rice singlet[20], which effectively supports a simpler one-band model 
with Cu only. Allais et ai.[23] proposed that the d-symmetry of these form factors, referred to as bond orders[21,22] because 
they are measured between the nearest neighbour Cu bonds, arise from the strong correlation but without other intertwined 
orders. Furthermore, there are also doubts regarding whether a strong correlation is present or even needed to understand 
of the superconducting mechanism[24]. However, the complexities of the phase diagram and some recent theoretical works 
have indicated the possibility of a new phase of matter, i.e., the pair density wave (PDW)[25-28], as discussed in detail in 
a recent review article[25]. The new states are considered to have intertwined orders of PDW and CDW or spin density 
waves(SDW)[25]. 



For quite some time, various calculations [29-39] on the Hubbard and t —J type models have revealed low-energy intertwined 
states appearing as stripes or bidirectional charge-ordered states, such as checkerboard(CB). However, these works usually 
involved different approximations and parameters, which often resulted in different types of charge-ordered patterns, and these 
studies were mostly concentrated at a hole concentration of 1 /8, which is the most notable concentration in early experiments. 
Hence, it is not clear if these results were the consequence of the invoked assumption or the approximation used, or if they are 
a generic results in the phase diagrams of cuprates. There were attempts to produce these CDWs or PDWs using a different 
approach, such as using a mean field theory to study a f — 7-like model but taking the strong correlation as only a renormalization 
effect of dispersion[21,22,40,41]. A spin-fluctuation mediated mechanism to produce these states was also proposed for the 
spin-fermion model[42]. Recently, a novel mechanism of PDW was proposed, i.e., Amperean pairing[28], by using the gauge 
theory formulation of the resonating-valence-bond picture. In most of these approaches, the wave vectors or periods of the 
density waves are related to special features of the Fermi surface, including nesting, hot spots or regions with large density of 
states. However, the opposite doping dependence of CDW periods, observed for 214 and \23(YBa2Cu^0^^s) compounds[10], 
makes the Fermi surface scenario worrisome. 

Amid all this confusion, recent numerical progress achieved by using the infinite projected entangled-pair states(iPEPS) 
method[43], has provided us with a new clue. It was found that the t — J model has several stripe states, with nearly degenerate 
energy as the uniform state and, with coexistent superconductivity and antiferromagnetism. When the number of variational 
parameters is extrapolated to infinity, the authors concluded that the anti-phase stripe, which has no net pairing, has slightly 
higher energy than the in-phase stripe with a net pairing, which in turn, also has slightly higher energy than the uniform 
state. This result is very consistent with the result of variational Monte Carlo calculations[29] based on the concept of the 
resonating-valence-bond picture[18]. Furthermore, the result is also consistent with that of renormalized mean-field theory by 
using a generalised Gutzwiller approximation(GWA)[44] to treat the projection operator in the t — J model[30,45]. Hence, the 
result provides strong support to more carefully examine the low energy states of the t — J model with the variational approach 
using GWA. 

Here, we report our findings from a much more extensive examination of the renormalized mean-field theory prediction 
using the GWA for the hole-doped t — J model. We find many unidirectional and bidirectional charge-ordered states with nearly 
degenerate energies as the uniform state, especially in the lightly doped regime; thus, it is a much more general phenomenon 
than previously thought. All of these states have intertwined orders of PDW, CDW and/or SDW. One of the CDW states, denoted 
as AP-CDW, reveals a bond order pattern with a much larger d-form factor than s' symmetry, as found in the experiment[16] 
with BSCCO {Bi2Sr2CaCu20%+x) and NaCCOC {Ca2-xNaxCu02Cl2)- Furthermore, just as in the experiment[17], it vanishes 
beyond 19% hole doping. However, not all these charge-ordered states have a dominant 7-form factor. For example, a different 
CDW intertwined with SDW and PDW, which is the familiar stripe reported long ago for 214 [1,2,29,31,36], has a larger s' 
form factor, as reported in the experiment[15]. We further show that this AP-CDW state could be easily altered to become 
a superconducting state with a global d-wave pairing symmetry, while locally, each bond does not have the perfect d-wave 
symmetry. Its spectra shows a large spatial variation at higher energies but with a d-wave nodal-like LDOS near zero energy as 
seen in the experiments [2,17]. 


Results and Discussions 

As mentioned above, the variational approach has been quite effective at capturing the physics of the strong correlation present 
in the t — J model. By using GWA, we can replace the strong constraint of forbidding the double occupancy of two holes on the 
same site in the variational wave function using Gutzwiller factors[32,33,44,45]. Then, one can use just mean field theory to 
find the various low energy states. Details about the calculation are discussed in the Methods section. 

In our mean field theory, there are four variational order parameters. Besides the hole density 5,, the local spin moment 
m] provides the antiferromagnetic correlation, the pair field Af ^ represents the local electron pairing order, and bond order 
Xjjfj is just the kinetic hopping term, where i is a site position and ij is the nearest neighbour bond. An iterative method is 
used to self-consistently solve the mean-field Hamiltonian //MfCEq. (S7) in the Supplementary Material(SM)) for all the 
parameters, of which there could be more than 60. The convergence is achieved for every order parameter if its value changes 
by less than 10^^ between successive iterations. All the calculations are performed on a 16 by 16 square lattice. To obtain 
various charge orders, specific patterns of 5,, mf, and Af ^ are input as initial values. The bond orders Xija always initially 
assumed to be uniform. In most cases, we will obtain only uniform solutions such as the d-wave superconducting(dSC) state 
and/or coexistent antiferromagnetic(dSC-AEM) state, but sometimes the states with charge-ordered patterns are found as a 
self-consistent solution. 
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Charge-ordered Patterns 

In addition to the two uniform solutions of a dSC state and a dSC-AFM state, there are many non-uniform charge-ordered 
states. For simplicity, we shall hrst present those charge-ordered states with a period of four lattice spaces (4ao), as listed in 
Table 1. Both the pair held and the spin moment mj’ could have positive and negative values. It turns out that if there is a 
SDW or a bidirectional spin CB (sCB) present, then it always has a period of with two domains of size 4ao with opposite 
antiferromagnetic directions joining together. The pair held has more choices. It could always be positive, with all of its x-bond 
pair held being positive and y-bond pair held being negative: thus, it would have a net total non-zero pair held. This is called an 
in-phase (IP) state, with a period of 4ao- However, just like the spin moment, the pair held could also have two domains with 
opposite signs and a domain wall in between: this state is known as the anti-phase (AP) state, with a period of 8ao- Thus, we 
could have four possible states for each unidirectional CDW or bidirectional charge CB (cCB), as we either have an IP or AP 
pair held with or without SDW. However, we only have three such states In Table 1 because we cannot hnd a solution with an 
IP pair held and CDW only. This result Is due to the choice of the commensurate period being 4ao- Later, we will show a state 
with a net pairing order or IP pairing state and CDW, which occurs if we do not require solutions to be commensurate with the 
lattice. 

Figure 1 shows a schematic illustration of the modulations of the pair held, charge density and spin moment for the three 
stripes with hole concentration of 0.1. The magnitude of the pair held is proportional to the width of the bond; red(cyan) 
denotes positive(negative) value. The size of the arrow is proportional to the spin moment, and the size of the circle represents 
the hole density. A similar hgure for the three bidirectional CB patterns is shown in Figure 7 in SM. There is one domain wall 
corresponding to the vanishing spin moment for IP-CDW-SDW in Figure la or the vanishing pair held for AP-CDW In Figure 
Ic. Both domain walls are present for the AP-CDW-SDW states in Figure lb. The hole density is always maximum at the 
domain wall with the vanishing spin moment. However, if there is no SDW, such as the AP-CDW stripe in Figure Ic, then the 
hole density is maximum at the domain wall with the vanishing pair held. This hndlng Is different from previous work without 
including the renormalized chemical potential[37]. 

Figure 2 shows energies as a function of hole concentration for all the states listed in Table 1 . The three unidirectional 
states are shown in the lower inset with blue triangles, circles, and diamonds representing IP-CDW-SDW, AP-CDW-SDW, and 
AP-CDW, respectively.The three CB states are shown in the upper inset with red triangles, circles and diamonds representing 
IP-cCB-sCB, AP-cCB-sCB, and AP-cCB, respectively. Unless specihcally mentioned, we only report site-centred results. 
Bond-centred solutions have essentially the same energies. The same results for the three CDW states were also reported in 
Ref. [30] at a 1/8 hole concentration. These mean-held GWA results are quite consistent with the numerical Monte Carlo result 
[29], which revealed that the uniform state has the lowest energy, followed by the in-phase stripe, and that the energy of the 
anti-phase stripe is slightly above that of both of them. However, the small energy differences are insignlhcant compared to the 
result of iPEPS[43], which showed the same ordering of states but with essentially degenerate energies. 

At approximately 12% doping in Figure 2, the spin moment becomes smaller, and the uniform dSC-AFM state merges 
into the dSC state. The difference from the original work of Ogata and Himeda[32,33], in which the spin moment vanished at 
10% doping, is due to the simplihed Gutzwiller factors used in Eq. (4). All the magnetic states, such as SDW and sCB, vanish 
at approximately 12% doping. The most surprising and important result shown in Eigure 2 is that in addition to the uniform 
dSC state, the AP-CDW state is most stable for a large doping range, from 0.08 to 0.18. The AP-cCB state also extends a 
little bit beyond the antiferromagnetic region. We only find the diagonal stripe state up to 6% doping. Another pattern that 
seems to be limited to small doping is IP-cCB-sCB, which is only found at doping less than 0.1. The general locations of 
these CB states in Eigure 2 are consistent with experimental observations that CB are seen more often at low doping [12,13]. 
Because the Gutzwiller factor gj in Eq. (4) is proportional to the hole density at the site, we expect the kinetic energy to be 
maximum at the domain wall (Eigure Ic), as shown in Table 2. Table 2 lists the values of hole density, the magnitude of the 
pairing order parameter and the kinetic energy K at each site, which are calculated by averaging the four nearest neighbour 
hopping amplitudes for AP-CDW at a 1/8 hole concentration. The kinetic energy and pairing order are calculated from the 
variational parameters x^ja respectively, by using Eq. (S9) in SM. Similar tables for other stripes and CB patterns are 

presented In Tables 3 and 4 in the SM. 

The red cross in Eigure 2 at the 1/8 hole concentration is the energy of a solution as a result of relaxing the requirement to 
have a commensurate 4ao period for the AP-CDW state. To alleviate the difficulty of considering incommensurate solutions in a 
hnite lattice calculation, we allow the state to have more than one single modulation period. In Eigure 3, the hole density, listed 
as the red numbers below the pattern, along with the magnitude of the pairing order parameter for both x and y bonds, listed in 
the top and bottom rows, are plotted along the direction of the modulation for a complex bond-centred stripe of length 16ao- 
It is very similar to the AP-CDW state. However, there is a remaining net constant d-wave pairing, with the system average 
Ax = —0.0056 and A^ = 0.0057. This mixture of the AP-CDW stripe with a small constant uniform pairing will produce a 
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d-wave nodal-like LDOS in addition to a PDW; hence, we have a nodal PDW or nPDW. There are several important results 
associated with the nPDW. Figure 3 shows that the hole density is indeed maximum at the domain walls near sites 4,7,10 and 
13. The maximum amplitude of pairing order A is about 0.03, which is roughly the same as adding the net pairing amplitude to 
that of the AP-CDW stripe in Table 2. It is most gratifying to observe that the d-wave pairing is globally maintained, although 
we have no way of controlling it during the iteration, with variables changing independently on each site. Contrary to the pure 
AP-CDW state without a net pairing, this state has a d-wave nodal spectrum at low energy, hence a nodal-like LDOS. In Figure 
4a, the LDOS of this stripe at 8 sites is plotted as a function of energy. The positions of these 8 sites are indicated in the inset of 
Figure 4a. The detailed LDOS at low energy is shown in Figure 4b. The large spatial variation of LDOS at high energies but 
always with a d-wave node near zero energy is quite consistent with the STM results in Ref. [3]. We have obtained this result by 
using a lattice of 16x16 supercells; please see the SM for details. 

A special feature of all these charge-ordered states is the large variation of the Gutzwiller factors from site to site. The 
values could change between nearest neighbours by a factor of 2 to 3. This unique property of strong correlated systems 

originates from the dependence on local hole density in the Gutzwiller factor, which is g\ = when we do not consider 

magnetic moments. This dependence on 5, is the consequence of being a Mott insulator when there are no doped holes. A slight 
variation of the hole density 5, will cause a large change in gj; in fact, dg\/d5i is proportional to g;/5,' ~ 1/v^. This factor 
dominates in the renormalized local chemical potential defined in Eq. (S6) when hole concentration is small. Thus, g\ is no 
longer a purely passive renormalization factor; now, it could alter the local chemical potential greatly and induce non-uniform 
charge orders. Although the factor associated with spin, in Eq. (4), is smaller, it also contributes to the local chemical 
potential. The strong susceptibility to the variation of local hole density makes a uniform state nnstable amidst inherent or 
extrinsic charge fluctuations. This effect is clearly more prominent in the lightly hole-doped regime, as demonstrated by the 
greater variety of charge-ordered states in the underdoped regime in Eigure 2. Another important effect of the Gutzwiller 
factor is that it introduces nonlinearity into the Bogoliubov-deGenne(BdG) equations(Eq. (S4)-(S6)), which can produce quite 
unexpected solutions. 

Bond Order 

So far, we have only discussed the pair held, hole density and spin moment; now, we shall consider more carefully the 
bond order Kij = \ Y.a{c\^cja) + (c^j^Cia)- The value of one-half in front of the summation is for averaging because there 
are two hopping terms for each bond. Now, it can be calculated by using the BdG solution and the Gntzwiller factor, i.e., 
Kij = 5 s'ijaXija s'jiaX'jia- Following the definition of bond order by Sachdev and collaborators [21,22,40] and Eujita et 
aJ.[16], by associating Kij+i ~ the tunneling current measured at the x-bond oxygen site can be obtained, similarly for 

the y-bond oxygen. The Fourier transform of these two quantities gives us the intra-unit-cell form factor. The Fourier transform 
of the AP-CDW state with a hole concentration of 1 /8 is schematically shown in Figure 5a. The size of the dot represents the 
magnitude; red (blue) represents a positive (negative) value. Because this is a 4ao stripe, in addition to values at Q = (0,0) and 
reciprocal lattice vectors denoted by the “+” sign, the modulation wave vector is (±7r/2ao,0), and the vectors are shifted by the 
reciprocal lattice vectors. The peaks at (±7r/2ao,0) are determined by Ay, while those at (±37r/2aojO) and (±7r/2ao,±27rao) 
are determined by Ajj- The ratio of Ao to Ay, or d/s', is approximately 7.5 in this case. This ratio is quite special for the 
AP-CDW state. For the IP-CDW-SDW stripe, the ratio is actually less than one. The schematic plots of the Fourier transform of 
IP-CDW-SDW and AP-CDW-SDW stripes are shown in Figure 8a and 8b in the SM, respectively. For the AP-CDW-SDW 
stripe, d/s' is approximately 1.2. The Fourier transform of the bond orders of the AP-cCB pattern is similar to that of AP-CDW 
with a dominant d-form factor, as discussed in the SM. 

The nPDW stripe shown in Figure 3 also has a large ^f-form factor with almost zero s'. The Fourier transform of its bond 
order is schematically shown in Figure 5b. The size of the dot scales with the magnitude of the c/-form factors, and red (blue) 
represents a positive (negative) value. The wave vector with a large amplitude is at 57r/8ao or its period is approximately 3.lag. 
This length is close to the separation between the domain walls of the pair held shown in Figure 3. The presence of smaller 
peaks at several wave vectors shows a mixture of different periods in the stripe. This result is expected if we add a constant 
pairing order to the AP-CDW stripe. 

Figure 5c is copied from Figure 3G of the STS work of Fujita et aJ.[16]. It shows the sum of real part of Fourier transform 
values of tunneling currents measured at Ox and Oy sites. Just like Figure 5a and Figure 5b, The value at (±3;r/2ao, 0) is larger 
than that at {±7:/2ao,0) and both have the same sign but opposite sign with respect to {±7:/2ao,3z2nao). In their sample there 
are two domains with density modulation in x and y directions, respectively. 

Another interesting result regarding the AP-CDW stripe is that its (7-form factor actually vanishes at an approximately 19% 
hole concentration, as shown in Figure 6 for both site-centred (blue dots) and bond-centred (red dots) solutions. We cannot 
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find the AP-CDW solution beyond 18% doping. This outcome is in excellent agreement with the results reported by Fujita 
et aJ.[17] in their Figure 3G which is copied as the inset of Figure 6. They measured the doping dependence of intensity of 
the modulation wave vector near (±3;r/2flo,0), which is associated with the density wave. The density wave disappeares at 
approximately 19% doping. Moreover, this 19% hole concentration is conspicuously close to the so-called quantum critical 
point[46]. We shall study this issue more in future work. 

Conclusion 

The results reported above are all based upon the well-established renormalized mean-field theory[45] and GWA[44] for a 
well-studied t — J model. Although they do not provide extremely accurate numbers, as many sophisticated numerical methods 
do, our results show that they do capture the most important physics of the strong correlation. This strong correlation provides 
a site-dependent Gutzwiller renormalization that produces many exotic solutions of PDW stripes and/or CBs intertwined 
with modulations of charge density and/or spin density. These results show quantitative agreement with some of the key 
experiments[3,12,13]. Because site-renormalization is extremely local, the effect of the Fermi surface or wave vectors kp is 
absent. Our model does not have the second or third neighbour hopping to provide a Fermi surface with nesting vectors or “hot 
spots” [21,40,46]. Thus, in our theory, there are no unique wave vectors for the charge density waves or CBs. Although we have 
mainly focused on the structures with a period of 4ao so far, our preliminary study also finds charge-ordered states with periods 
of 5ao and even 3ao. States with a longer period should be possible, and they could also have degenerate energies[34,43]. If we 
allow a pattern with multiple periods, such as the nPDW stripe shown in Figures 3 and 5b, we could have states with fractional 
or incommensurate periods. A detail study of all these will be conducted in the future, as well as a study of the effect of having 
values of J/t away from 0.3. 

An important conseqnence of having all these charge-ordered states originating from the same Hamiltonian and physics 
is that these states are not the usual “competing states” we are familiar with. They do not stay in a deep local minima in the 
energy landscape. They are actually quite fragile and can easily evolve into each other, as we have already demonstrated with 
the nPDW stripe, which evolved from a mixture of AP-CDW and an uniform d-SC state. Other examples of the mixmre of 
stripes listed in Table 1 can be easily constmcted. For real cuprates, there are many other interactions in addition to our t and J 
that will alter the preferences of these states. For example, a weak electron lattice interaction could make the IP-CDW-SDW 
stripe much more stable against the dSC-AFM state[36]. Including special Fermi surface features could also enhance CDW for 
certain periods. However, none of these interactions are as important and necessary as the site renormalization due to strong 
Mott physics to produce these charge-ordered states. The effect of finite temperature will certainly bring in the entanglement 
of these states and much more complicated phenomena, snch as psendogap. Developing a method for generalising GWA to 
include the temperature effect remains as a big challenge. 


Methods 


We introduce the t — J Hamiltonian[18] on a square lattice of Cu by using 

H = - Y. PGticlcja+H.C.)PG + Y JSi ■ Sj 

{U),o {ij) 

where nearest neighbour hopping f, as our energy unit, is set to 1, and J is set to 0.3. Fb = Hill ~ *^he Gutzwiller 

projection operator, while n,cj = cj^c/a stands for the number operator for site i. Spin a is equal to ±. 5, is the spin one-half 
operator at site i. The Fermi surface of the uniform state is quite simple, without nesting parts, and does not intersect with the 
magnetic Brillouin zone boundary, thus avoiding hot spots. 


Following the idea of Gutzwiller[44] and work of Himeda and Ogata[32,33], we replace the projection operator(P( 3 ) with 
the Gutzwiller renormalization factors. The renormalized Hamiltonian now becomes 


^ = - E g'jat(clcja+H.C.) +Y-^ 

(iJ) 




( 2 ) 


where and gfj^ are the Gutzwiller factors, which are dependent on the values of local AF moment mj’, pair field 

bond order Xija’ density 5,; 


mj* = (‘Fo|5j|‘Fo) 

Al, = (j('Fo|c,aC;s|‘Po) 

Xija = ('PokLc.aln) 
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where I'Pq) is the unprqjected wavefunction. The superscript v is used to denote that these quantities are different from the real 
physical quantities for comparison with the experiments. Their relationship is given in Eq. (S9). As for the Gutzwiller factors, 
we follow the work of Yang et al.[30]; they used a slightly simplified version of Ogata and Himeda[32,33], which was also used 
by Christensen et aJ.[34]. The factors are given as 


Si jo SioSjo 


Sia 


25,(1-5,) l + 5i + a2mj 


1 — 5? + 1 + 5,' — (72mJ’ 


^,xy 5,rv 

SiJ =Si 8j 

2(1-5,) 


s.xy 

8i = 


1 — 5? +4(mj’)2 


s,z s.xy 
8ii =8i-- 


Xij = l- 




''ij 


I jij 


2((A:))2 + (zf)2)-4m>;; 
12(1-5,■)(l-5,)((AV.)2 + (i,f)2) 


1^(1 - 5,2 +4(ot;') 2)(1 - 5? +4(mp2) 


(4) 


where Af- = Y.a^'ija/'^ ^tid xjj = T,oXijal^- the presence of antiferromagnetism, A'^j^ ^ AV.^. The derivation of the 
mean-field self-consistent equations is described in the SM. 
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Figure Legends 

Figure 1 

Schematic illustration of modulations for stripe like patterns: (a)lP-CDW-SDW (b)AP-CDW-SDW (c)AP-CDW respectively. 
Size of the circle represents the hole density. The width of the bond around each site represents the amplitude of pairing 
A(A L<t A<t) sign is positive(negative) for red(cyan). The size of black arrows represents the spin moment. The average 

hole density is about 0.1. 

Figure 2 

Energy per site as a function of hole concentration. Six states are shown in the main figure with notations defined in Table 2. The 
lower(upper) inset is for stripe(CB) patterns. Blue triangles, circles, and diamonds are for IP-CDW-SDW, AP-CDW-SDW, and 
AP-CDW respectively. And red triangles, circles and diamonds are for IP-cCB-sCB, AP-cCB-sCB, and AP-cCB respectively. 

Figure 3 

Schematic illustration of modulations for nPDW stripe. The numbers in red denote the hole dnesity at each site while the 
numbers in black below them represent the pairing amplitude in y direction. The rest numbers above the figure stand for the 
pairing amplitude in x direction. Here our pairing amplitudes denote {{cj^Cji)). Note that in this figure neither the size of 
circles nor the width of bonds represent amplitudes. The hole concentration is 0.125. 

Figure 4 

(a) LDOS at 8 sites plotted from energy 0.6t to -0.6t. The inset shows hole density along the modulation direction of the nPDW 
stripe and (b) from 0.2t to -0.2t but shifted vertically for clarity. 

Figure 5 

Schematic illustration of the Fourier transform of bond orders of (a) AP-CDW state and (b) the nPDW stripe in a lattice of 
16ao*16ao- ”+” signs are at the four reciprocal lattice vectors (±27r/aO)0) and (0,±27r/ao) and their nearby medium size dots 
are shifted from them by (±7r/2ao,0). The large dot at center isQ— (0,0) and has two red small dots nearby at (±7r/2aojO). 
The inner dotted square is the boundary of first Brillouin zone, (c) is copied from Figure 3G of the STS work of Fujita et aJ.[16]. 
It shows the sum of real part of Fourier transform values of tunneling currents measured at Ox and Oy sites. Unlike (a) and (b) 
that only has one domain of density modulation in the x direction, this sample has two domains with both x and y direction 
modulations. 

Figure 6 

Magnitude of the d-form factor for the AP-CDW stripe as a function of doped hole concentration. Blue dots are for site-centered 
AP-CDW stripe and red ones for bond-centered AP-CDW. The inset is copied from Figure 3G of the STS work of Fujita et 
aJ.[17] showing the doping dependence of intensity of the modulation wave vector near (±3;r/2flo,0), which is associated with 
the density wave. This modulation vanishes at 19% hole concentration. 
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pair field 

charge modulation 

spin modulation 

IP-CDW-SDW 

in-phase 

stripe 

yes 

AP-CDW-SDW 

anti-phase 

stripe 

yes 

AP-CDW 

anti-phase 

stripe 

zero 

IP-cCB-sCB 

in-phase 

checkerboard 

yes 

AP-cCB-sCB 

anti-phase 

checkerboard 

yes 

AP-cCB 

anti-phase 

checkerboard 

zero 

dSC 

uniform 

uniform 

zero 

dSC-AFM 

uniform 

uniform 

uniform 

diag 

in-phase 

stripe along (1,1) 

yes 


Table 1 . Definition of various nearly degenerate states with respect to the intertwined orders: pair field, charge density, and 
spin moment. Besides the two uniform solutions, d-wave superconducting (dSC) state and coexistent antiferromagnetic 
(dSC-AFM) state, all the states to be considered in this paper, unless specifically mentioned, have modulation period 4ao for 
charge density and bond order. IP(AP) means the pair field is in-phase with period 4ao (anti-phase with period 8ao)- IP has a 
net pairing order and AP has none. SDW is the spin density wave with period Saq. sCB (cCB) denotes the checkerboard pattern 
of spin (charge) and diag means the diagonal stripe which has in-phase pair field and spin modulation. 


site number 

1 

2 

3 

4 

5i 

0.1315 

0.1256 

0.1168 

0.1256 

A; 

0 

0.0194 

0.0247 

0.0194 

K, 

0.092 

0.0866 

0.0799 

0.0866 

Kij+f 

0.1151 

0.0901 

0.0625 

0.0901 

KiJ+x 

0.0688 

0.0972 

0.0972 

0.0688 


Table 2. Hole density and order parameters at each site for an AP-CDW stripe at 0.125 doping. A, is the average of pairing 
order of the four bonds at site i. Ki is the average kinetic energy at each site and (f^;,!+.v) are the bond orders in the y (x) 
direction. These parameters are calculated according to Eq. (S9). 
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Following the renormalized mean-field theory[l] by using the GWA as in the works of Yang et a].[2], we derive the formula 
we used for solving the BdG equations. After we replace the projection operator by the Gutzwiller factors and use the mean-field 
order parameters defined in Eq.(3), the energy of the renormalized Hamiltonian(Eq.(2)) becomes 


£ = ('Poii/i'Po)= - 

ij-(y 


+ //,c,)- Ey(f 




A"*- 

°‘j ij<y 


2 A' 


ij(j 


A” A*"- 






8ii Z, 


‘j(^ 




2 Z 


‘j<^ 


(i-j) 


(SI) 


Next we want to minimize the energy under two constraints:^,-n, = Ne and (‘PqI'Pq) = 1. Thus our target function to be 
minimized is 


IT = ('Po|//|'Po) - A((‘Po|‘Po) -l)-^{Y^ni-Ne) 

i 

The mean-field Hamiltonian now becomes 

HmF - 2^ 3 „ C^^Cja+H.C. -f 2^ aCiaCja+H.C.+^-^—tlia 

i.j.a°Xija i,a 


(S2) 


(S3) 
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Eq.(S3) satisfies the Schrodinger equation = ■^I'Eo)- The three derivatives are defined as 
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and the effective local chemical potential is defined as 
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y. dW dg^f 

id^dma 


dW dg\ 


ija' 




(S4) 


(S5) 


(S6) 


The last term dg’.j^,/dni(j in the effective local chemical potential gives the biggest contribution. This was not included in 
previous works using GWA to study CB[3,4], and their results have different patterns compared with ours. In addition the 
energy variation between different charge-ordered states is much larger than our nearly degenerate results. 

Now Hmf can be rewritten as BdG equations, 



We can diagonalize the Hmf to obtain equal number of positive and negative eigenvalues with their corresponding eigenvectors 
{u",\X). With these eigenvectors, we can determine the order parameters at zero temperature by following equations 

ri- 

”4 = (c4Gi)o = El''^l^ 

n+ 

"V (S8) 

A4 = -(c4Cit)o = E«>r 

H + 

X 4 t = ( 4 ";t)o = E“>r 

ri- 

xhi={c]ic,i)«='L''yr 


The sum for n+(n_) means the set of eigenvectors with positive(negative) energies. An iterative method is used to solve Hmf 
self-consistently. The convergence is achieved for every order parameter if its value changes less than 10^^ between successive 
iterations. After the self-consistency is achieved, we calculate order parameters, their formula are 


A; E^^f <r^i+4<rAi,r+f,(7 T g; (j'?;—.fjcA; ,'_;f (j 8i,a8i+y,a^i,i+y,a 8i,a8i—y,(y^i,i—y,a)l^'’ 


=( ^+V 

^ (7 

^i,i+y ^j8i,i+y,a {^ia^i+yo ) T 8i+y^i^a {^i+ya^i<x ) i 
^ a 

Ki = {Kij+y + Ki i^y -\- + /T;/_ji) /4 


(S9) 


The values for the above quantities was shown in Table 2 of the main text for a typical AP-CDW stripe. Here we show the 
values for two other stripes in Table 3 and three CB patterns in Table 4. 

A schematic illustration of CB like patterns is shown in Figure 7. Definitions of symbols are same as Figure 1 in the main 
text. Again, same as stripes shown in Figure 1 from the main text, we have the maximum hole density at sites either on the 
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domain walls of AFM or pair field if AFM is absent. The latter part is different from previous results using GWA to study 
CB[3,4]. 

We can also examine the symmetry of bond orders, and Ki i+^, as we did in Figure 5 of the main text by examining 
the Fourier transform. CB in Figure 9 shows clearly it can be thought of as the linear combination of stripes in x and y direction. 
The small dots inside the dotted square or the first Brillouin zone is proportional to the s' form factor or Ay discussed in the 
main text, and the outside dots are related to the d-form factor or Aj). Just like AP-CDW stripe, the AP-cCB in Figure 9c also 
has a much larger ratio of d/s'. For IP-CDW-SDW stripe and IP-cCB-sCB, the ratio is less than one. For AP-CDW-SDW and 
AP-cCB-sCB, the ratio is about 1. 

We can also use the BdG solutions to calculate LDOS as shown in Figure 4 of the main text for the nPDW stripe. Here 
we use the supercell method[5] to calculate LDOS. Each cell has Nx x Ny lattice points and we have = Mx x My cells. We 
can now reduce the Hamiltonian from 2MxNx x IMyNy to Mx x My matrix equations each with dimension 2Nx x 2Ny. LDOS is 
calculated by the equation 

P'(^)=i^E[(^L)'l<(K)l"5(£-£„(K)) + (4)2|v«(K)|25(£+£„(K))] (SlQ) 

K,„ 

where K = 2n{j^^, nx G [OjMj: — 1] and ny G [0,My — 1]. Also we replace the delta function by a Lorentzian function 

with the width set to be O.Olt in this paper. In Figure 4 of the main text, we had used Mx = My =Nx=Ny= 16. 
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IP-CDW-SDW 




AP-CDW-SC 

»W 

site number 

1 

2 

3 

4 

site number 

1 

2 

3 

4 


0.1141 

0.1252 

0.1369 

0.1252 


0.1071 

0.1268 

0.1396 

0.1268 

l»u| 

O.lCl 

0.1146 

0 

0.1146 

l»ul 

0.2315 

0.1189 

0 

0.1189 

A. 

0.0255 

0.0256 

0.0256 

0.0256 

A. 

0 

0.0219 

0.0273 

0.0219 

A'. 

0.0771 

0.0844 

0.0925 

0.0844 

A'. 

0.0726 

0.0844 

0.0926 

0.0844 


0.0723 

0.0856 

0.1003 

0.0856 


0.0971 

0.0986 

0.0927 

0.0986 


0.0818 

0.0846 

0.0846 

0.0818 


0.048 

0.0924 

0.0924 

0.048 


Table 3. Values of several order parameters for IP-CDW-SDW and AP-CDW-SDW stripes at 0.125 doping. 
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IP-cCB-sCB 



1 

2 

3 

4 

1 

0.0318 

0.0594 

0.1037 

0.0594 

2 

0.0.594 

0.0849 

0.1205 

0.0849 

3 

0.1037 

0.1205 

0.1378 

0.1205 

4 

0.0594 

0.0849 

0.1205 

0.0849 


Ai 

1 

2 

3 

4 

1 

0.0164 

0.0162 

0.0157 

0.0162 

2 

0.0162 

0.0167 

0.0172 

0.0167 

3 

0.0157 

0.0172 

0.0189 

0.0172 

4 

0.0162 

0.0167 

0.0172 

0.0167 


\mi\ 

1 

2 

3 

4 

1 

0.3614 

0.2879 

0 

0.2879 

2 

0.2879 

0.2271 

0 

0.2271 

3 

0 

0 

0 

0 

4 

0.2879 

0.2271 

0 

0.2271 


AP-cCB-sCB 


Ki 

1 

2 

3 

4 

1 

0.0216 

0.0403 

0.0683 

0.0403 

2 

0.0403 

0.0567 

0.0793 

0.0.567 

3 

0.0683 

0.0793 

0.0933 

0.0793 

4 

0.0403 

0.0567 

0.0793 

0.0.567 


A; 

1 

2 

3 

4 

1 

0 

0 

0 

0 

2 

0 

0.0051 

0.0042 

0.0051 

3 

0 

0.0042 

0.0055 

0.0042 

4 

0 

0.0051 

0.0042 

0.0051 


Si 

1 

2 

3 

4 

1 

0.0686 

0.083 

0.1106 

0.083 

2 

0.083 

0.0959 

0.1176 

0.0959 

3 

0.1106 

0.1176 

0.1257 

0.1176 

4 

0.083 

0.0959 

0.1176 

0.0959 


|m,| 

1 

2 

3 

4 

1 

0.3793 

0.2857 

0 

0.2857 

2 

0.2857 

0.221 

0 

0.221 

3 

0 

0 

0 

0 

4 

0.2857 

0.221 

0 

0.221 


AP-cCB 


lu 

1 

2 

3 

4 

1 

0.0464 

0.0573 

0.0757 

0.0573 

2 

0.0573 

0.0671 

0.0819 

0.0671 

3 

0.0757 

0.0819 

0.0906 

0.0819 

4 

0.0573 

0.0671 

0.0819 

0.0671 


Si 

1 

2 

3 

4 

1 

0.1229 

0.1109 

0.0926 

0.1109 

2 

0.1109 

0.1039 

0.0884 

0.1039 

3 

0.0926 

0.0884 

0.078 

0.0884 

4 

0.1109 

0.1039 

0.0884 

0.1039 


Ki 

1 

2 

3 

4 

1 

0.0865 

0.0786 

0.0649 

0.0786 

2 

0.0786 

0.0717 

0.0601 

0.0717 

3 

0.0649 

0.0601 

0.0518 

0.0601 

4 

0.0786 

0.0717 

0.0601 

0.0717 


A, 

1 

2 

3 

4 

1 

0 

0 

0 

0 

2 

0 

0.01.54 

0.0198 

0.0154 

3 

0 

0.0198 

0.0252 

0.0198 

4 

0 

0.01.54 

0.0198 

0.0154 


Table 4. Values of several order parameters for checkerboard patterns. The hole concentration is 0.1 for AP-cCB-sCB and 
AP-cCB but 0.09 for IP-cCB-sCB. 
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Figure 7. Schematic illustration of modulations for CB like patterns: (a)IP-cCB-sCB (b)AP-cCB-sCB (c)AP-cCB 
respectively. Definitions of all symbols are same as Figure 1 of the main text. The average hole density is 0.1 for AP-cCB-sCB 
and AP-cCB and 0.09 for IP-cCB-sCB. 


17/19 

































Figure 8. Schematic illustration of the Fourier transform of the bond orders for (a)IP-CDW-SDW stripe and for 
(b)AP-CDW-SDW. The dot size scales with the magnitude and red (blue) for positive(negative) values. ”+” signs are at the four 
reciprocal lattice vectors (±2;r/ao,0) and (0,±2;r/ao) and their nearby medium size dots are shifted from them by 
(±7r/2aojO). The center large dot '\sQ— (0,0) and has two red small dots nearby at (±7r/2ao,0). The inner dotted square is 
the boundary of first Brillouin zone. The doping for both stripes is 1/8. 
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Figure 9. Schematic illustration of the Fourier transform of the bond orders for CB patterns (a)IP-cCB-sCB, (b)AP-cCB-sCB, 
and (c)AP-cCB. The hole density is (a)5 = 0.09, (b)5 = 0.1, and (c)5 =0.1. All the dots are shifted from Q — (0,0) and the 
four reciprocal lattice vectors (denoted by ”+” sign) by (±2;r/ao,0) or (0,±2;r/flo)- The notations are the same as those in 
Figure 8. 
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